#coding:utf-8

import random
import numpy as np
import pylab as p

class ListaPontosNM:
	def __init__(self,total):
		self.lista_pt = [0]*total
		
		self.total = total
		self.pt_medio = 0
		
	def insert_ponto(self,pos,pt):
		self.lista_pt[pos] = (pt,F(pt))
		
	def replace_ponto(self,pos,pt):
		self.insert_ponto(pos,pt)
		
	def reduz_simplex(self):
		for i in range(1,self.total):
			self.replace_ponto(i, (self[i][0]+self[0][0])/2 )		

	def __getitem__(self,i):
		return (self.lista_pt[i])
		
	def sort(self):
		self.lista_pt.sort(key=lambda lpnm:lpnm[1])
		
		self.pt_medio = np.zeros( len(self.lista_pt[0][0]) )
		for i in range(1,self.total-1):
			self.pt_medio+= self.lista_pt[i][0]
			
		self.pt_medio = self.pt_medio/(self.total-1)
		

def teste(a,b,c,x,y,z,r1,r2):
	return pow((x-a),2)+pow((y-b),2)+pow((z-c),2)-pow((r1+r2),2)

def gen_pontos(nvar):
	zeros = np.zeros( (nvar+1,nvar) )
	pts = ListaPontosNM(nvar+1)
	for i in range(0,nvar):
		zeros[i,i] = 15
		pts.insert_ponto(i,zeros[i])
		
	for i in range(0,n):
		ind = 3*i
		zeros[nvar,ind] = 30
		zeros[nvar,ind+1] = 30
		zeros[nvar,ind+2] = 60
		
	pts.insert_ponto(nvar, zeros[nvar])
	
	return pts
	
def gen_esferas(n):
	xi = [ 100*random.random() for j in range(0,n) ]
	yi = [ 100*random.random() for j in range(0,n) ]
	zi = [ 100*random.random() for j in range(0,n) ]
	
	ri = [ int( 100*random.random() ) for j in range(0,n) ]
	
	return (xi,yi,zi,ri)
	
class Util:
	MAX_ITER = 500
	FATOR_1 = 0
	FATOR_2 = 1

n=5

xi = [2,2,2]
yi = [2,2,2]
zi = [2,2,2]
ri = [6, 7, 2, 3, 4]

'''	
res = gen_esferas(n)
xi = res[0]
yi = res[1]
zi = res[2]
ri = res[3]
'''


def soma_v(v1,v2):
	res = [0]*len(v1)
	for i in range(0,len(v1)):
		res[i] = v1[i]+v2[i]
	return res
	
def interno_v(v1,v2):
	res = 0
	for i in range(0,len(v1)):
		res += v1[i]*v2[i]
		
	return res	
	
def escalar_v(v1,k):
	res = [0]*len(v1)
	for i in range(0,len(v1)):
		res[i] = v1[i]*k
		
	return res	
	
def norma(v):
	res = 0
	for i in range(0,len(v)):
		res+=pow(v[i],2)
	res = np.sqrt(res)
	return res

	
def Ei(i,v):
	ind_i=3*i
	res = 0
	for k in range(0,n):
		ind_k=3*k			
		
		if i==k:
			continue
			
		parc = pow( (v[ind_i]-v[ind_k]),2 ) + pow( (v[ind_i+1]-v[ind_k+1]),2 ) + pow( (v[ind_i+2]-v[ind_k+2]),2 ) - pow( (ri[i]+ri[k]),2 )
		res+=pow(min(parc,0),2)
		
	return res
	
def Fi(i,v):
	ind =3*i
	res = Util.FATOR_1*pow(v[ind],2)+pow(v[ind+1],2)+pow(v[ind+2],2)
	res += Util.FATOR_2*Ei(i,v)
	
	return res
	
def F(v):
	res = 0
	
	max = [-100,-100,-100]
	min = [100,100,100]	
	for i in range(0,n):
		res+=Fi(i,v)
		
		ind = 3*i
		if v[ind]+ri[i] > max[0]:
			max[0] = v[ind]+ri[i]
		if v[ind+1]+ri[i] > max[1]:
			max[1] = v[ind+1]+ri[i]
		if v[ind+2]+ri[i] > max[2]:
			max[2] = v[ind+2]+ri[i]			
			
		if v[ind]-ri[i] < min[0]:
			min[0] = v[ind]-ri[i]
		if v[ind+1]-ri[i] < min[1]:
			min[1] = v[ind+1]-ri[i]
		if v[ind+2]-ri[i] < min[2]:
			min[2] = v[ind+2]-ri[i]
		
	res += ( ( abs(max[0]-min[0]) )*( abs(max[1]-min[1]) )*( abs(max[2]-min[2]) ) )
		
	return res

def volume(v):
	max = [-100,-100,-100]
	min = [100,100,100]	
	
	for i in range(0,n):
		ind=3*i
		if v[ind]+ri[i] > max[0]:
			max[0] = v[ind]+ri[i]
		if v[ind+1]+ri[i] > max[1]:
			max[1] = v[ind+1]+ri[i]
		if v[ind+2]+ri[i] > max[2]:
			max[2] = v[ind+2]+ri[i]			
			
		if v[ind]-ri[i] < min[0]:
			min[0] = v[ind]-ri[i]
		if v[ind+1]-ri[i] < min[1]:
			min[1] = v[ind+1]-ri[i]
		if v[ind+2]-ri[i] < min[2]:
			min[2] = v[ind+2]-ri[i]
		
	return ( ( ( abs(max[0]-min[0]) )*( abs(max[1]-min[1]) )*( abs(max[2]-min[2]) ) ),min,max )
	
	
def main():
	epson = 1e-4
	it = 0
	
	nvar = 3*n
	pts = gen_pontos(nvar)
	parada = 100
	
	while(parada>epson and it<Util.MAX_ITER):
		print it,parada
		pts.sort()
		xr = pts.pt_medio + (pts.pt_medio-pts[nvar][0])
		f_xr = F(xr)
		
		if f_xr < pts[0][1]:
			xe = pts.pt_medio + 2*(pts.pt_medio-pts[nvar][0])		
			f_xe = F(xe)
			
			
			if f_xe < f_xr:
				#print 11
				pts.replace_ponto(nvar,xe)
			else:
				#print 12
				pts.replace_ponto(nvar,xr)
		
		else:
			for i in range(1,nvar):
				if f_xr < pts[i][1]:	
					#print 21
					pts.replace_ponto(i,xr)
					break
			else:
				if f_xr < pts[nvar][1]:
					xce = pts.pt_medio + 0.5*(pts.pt_medio-pts[nvar][0])
					f_ce = F(xce)
					
					if f_ce < f_xr:
						#print 31
						pts.replace_ponto(nvar,xce)
					else:
						#print 32
						pts.replace_ponto(nvar,xr)
				else:					
					xci = pts.pt_medio - 0.5*(pts.pt_medio-pts[nvar][0])
					f_ci = F(xci)
					
					if f_ci < pts[nvar][1]:
						#print 41
						pts.replace_ponto(nvar,xci)
					else:
						#print 42
						pts.reduz_simplex()
		
		parada = 0
		f_xm = F(pts.pt_medio)
		
		
		for i in range(0,pts.total):
			parada+= pow( norma(pts[i][0] - pts.pt_medio),2 )
		parada/=pts.total
		parada = np.sqrt(parada)
		
		
		'''
		for i in range(0,pts.total):
			parada+= pow( (f_xm-pts[i][1]),2 )
		parada/=pts.total
		'''
		
		it+=1
	
	
	print it
	return pts.pt_medio
	
	
x = main()


#Print Scilab plot

print "v=[",
for i in range(0,len(x)):
	if i%3==0 and i!=0:
		print ";",
	print x[i],
print "]"

print "ri=",ri

resp = volume(x)
print "caixa=[ ",resp[1],"],[",resp[2],"]"
print "Volume=",resp[0]
